home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_10.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.2 KB  |  147 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.10. (Finite-Difference Method).
  9. % Section    9.9,    Finite-Difference Method, Page 496
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the finite difference
  13. % method for solving the boundary value problem
  14.  
  15. %       x`` = p(t)x`(t) + q(t)x(t) + r(t)    
  16. % with  x(a) = alpha, x(b) = beta
  17.  
  18. %    Store p(t),q(t), r(t)    in the M-files p.m  q.m  r.m 
  19. % function z = p(t)
  20. % z = 2*t/(1 + t^2);
  21.  
  22. % function z = q(t)
  23. % z = - 2/(1 + t^2);
  24.  
  25. % function z = r(t)
  26. % z = 1;
  27.  
  28. delete p.m
  29. diary p.m; disp('function z = p(t)');...
  30.            disp('z = 2*t/(1 + t^2);');...
  31. diary off;
  32.  
  33. delete q.m
  34. diary q.m; disp('function z = q(t)');...
  35.            disp('z = - 2/(1 + t^2);');...
  36. diary off;
  37.  
  38. delete r.m
  39. diary r.m; disp('function z = r(t)');...
  40.            disp('z = 1;');...
  41. diary off;
  42.  
  43. % Remark. p.m  q.m  r.m findiff.m trisys.m are used for Algorithm 9.10
  44. p(0); q(0); r(0);
  45. pause    % Press any key to continue.
  46.  
  47. clc;
  48. %    Place the endpoints of [a,b] in  a  and  b.
  49.  
  50. %    Place the initial  value x(a) in  alpha.
  51.  
  52. %    Place the terminal value x(b) in  beta.
  53.  
  54. %    Place the number of subintervals in  n.
  55.  
  56. a  = 0;
  57. b  = 4;
  58. alpha =  1.25;
  59. beta  = -0.95;
  60. n  = 20;
  61.  
  62. [T,X] = findiff('p','q','r',a,b,alpha,beta,n);
  63. points = [T;X];
  64.  
  65. pause    % Press any key to see the solution points.
  66.  
  67. clc; clg;
  68. c = min(X)-0.2;
  69. d = max(X)+0.2;
  70. axis([a b c d]);...
  71. plot(T,X,'-g');...
  72. hold on;...
  73. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  74. if n<=30,
  75.   plot(T,X,'or');
  76. end;...
  77. xlabel('t');...
  78. ylabel('x');...
  79. Mx1 = 'The finite difference solution to x`` = f(t,x,x`).';...
  80. title(Mx1);...
  81. grid;...
  82. axis;...
  83. hold off;...
  84. shg; pause    % Press any key to continue.
  85.  
  86. Mx2 = '     t(k)               x(k)'; clc;
  87. clc,echo off,diary output,...
  88. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  89. diary off,echo on
  90.  
  91. pause    % Press any key to perform extrapolation.
  92.  
  93. %    Place the endpoints of [a,b] in  a  and  b.
  94.  
  95. %    Place the initial  value x(a) in  alpha.
  96.  
  97. %    Place the terminal value x(b) in  beta.
  98.  
  99. %    Place the number of subintervals in  n.
  100.  
  101. a  = 0;
  102. b  = 4;
  103. alpha =  1.25;
  104. beta  = -0.95;
  105. n  = 20;
  106.  
  107. [T,X1] = findiff('p','q','r',a,b,alpha,beta,n);
  108. [T,X2] = findiff('p','q','r',a,b,alpha,beta,2*n);
  109. X2 = X2(1:2:length(X2));
  110. [T,X3] = findiff('p','q','r',a,b,alpha,beta,4*n);
  111. X3 = X3(1:4:length(X3));
  112. T = T(1:4:length(T));
  113. Z1 = (4*X2-X1)/3;
  114. Z2 = (4*X3-X2)/3;
  115. X = (16*Z2-Z1)/15;
  116. points = [T;Z1;Z2;X];
  117.  
  118. pause    % Press any key to see the solution points.
  119.  
  120. clc; clg;
  121. c = min(X)-0.2;
  122. d = max(X)+0.2;
  123. axis([a b c d]);...
  124. plot(T,X,'-g');...
  125. hold on;...
  126. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  127. if n<=30,
  128.   plot(T,X,'or');
  129. end;...
  130. xlabel('t');...
  131. ylabel('x');...
  132. Mx1 = 'The finite difference solution to x`` = f(t,x,x`).';...
  133. title(Mx1);...
  134. grid;...
  135. axis;...
  136. hold off;...
  137. shg; pause    % Press any key to continue.
  138.  
  139. format long
  140. Mx2 = 'Extrapolated solutions z1(k), z2(k) and x(k).';
  141. Mx3 = '     t(k)               z1(k)              z2(k)              x(k)';
  142. clc,echo off,diary output,...
  143. disp(''),
  144. disp(''),disp(Mx1),disp(''),disp(Mx2),...
  145. disp(''),disp(Mx3),disp(points'),...
  146. diary off,echo on
  147.